function covm = covmTS(R, HL, annual_factor)
% covmTS: Estimate time varying covariance matrix
% covm = covmTS(R, halflife, annual_factor)  Estimate time varying covariance matrix
% using simple expanding window if no halflife provided.  EWMA otherwise.
% Example:
%    [V] = covmTS(R, halflife, annual_factor);
%
% TODO: add separate vol and corr halflifes.  demeaning in a sensible way.
% consider capping outliers.  Turn daily return data into monthly covm TS 
% structure rather than daily covm TS.
%

N = length(R.dates);


% Build initial covm using halflife
if exist('HL','var')
    gamma = 0.5^(1/HL);
    decaypwr = (HL-1:-1:0)';
    wgts = (gamma).^decaypwr;
    wgts = wgts/sum(wgts);
    init_covm = R.data(1:HL,:)'*diag(wgts)*R.data(1:HL,:);
else
    error('no functionality for missing halflife yet.');
end

% Build covariance matrix time series
covm = buildTS([], R.header, R.dates, '3D');  % Covariance matrix

last_covm_t = init_covm;
for t = 1 : N

   rt = R.data(t, :)'; % Single period returns

   % Estimate covariances through time
   covm_t = gamma * last_covm_t + (1 - gamma) * (rt * rt');

   % Update t-1 values
   last_covm_t = covm_t;
   
   % Save current time period outputs
   covm.data(:, :, t) = covm_t;
   
end

% Annualize
if exist('annual_factor','var')
    covm = multTS(covm, annual_factor);
end
% Lag cov matrix for use in backtesting
covm.data = cat(3, covm.data(:, :, 1), covm.data(:, :, 1:end-1));
covm.data = cat(3, covm.data(:, :, 1), covm.data(:, :, 1:end-1));

return
